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ABSTRACT 

The r-decay library TAUOLA has gained popularity over the last decade. However, with 
the continuously increasing precision of the data, some of its functionality has become 
insufficient. One of the requirements is the implementation of decays into five mesons 
plus a neutrino with a realistic decay amplitude. This note describes a step into this 
direction. For the 27r~7r + 27r° mode the three decay chains r~ — > a{v — > p~[-^ vr _ 7r°)a;(^ 
7r~7r + 7T°)z/, t~ — > a\v — > a 1 "(— > 27r~7r + )/o(^ 2ir°)v, and r~ — »• ajfv — > ajf(— > 7r~27r°)/o( — *■ 
7T + 7r— )v are introduced with simple assumptions about the couplings and propagators of 
the various resonances. Similar amplitudes (without the pu contributions) are adopted 
for the 7r~47r° and 37r~27r + modes. 

The five-pion amplitude is thus based on a simple model, which, however, can be 
considered as a first realistic example. Phase-space generation includes the possibility of 
presampling the u and a\ resonances, in one channel only, however. This is probably 
sufficient for the time being, both for physics applications and for tests. 

The technical test of the new part of the generator is performed by comparing Monte 
Carlo and analytical results. To this end a non-realistic, but easy to calculate, purely 
scalar amplitude for the decay into five massless pions was used. 
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1. INTRODUCTION 



Early studies of semileptonic tau decays have concentrated on final states with few 
mesons only. The recent advent of high-statistic samples in experiments at LEP and 
CESR, and the perspective of further increasing event rates at B-meson factories allow and 
require the study of relatively rare decay modes and thus of multibody final states. These 
configurations are of particular importance for a precise determination of the semileptonic 
branching ratio and for an improved limit on the mass of the tau neutrino. Final states 
with up to six pions have been observed up to now. 

Complementary to the experimental studies, Monte Carlo simulations are required 
to determine the efficiencies of the detectors. The generator TAUOLA [1-3] has been 
specifically designed to simulate a wide variety of tau-decay modes and includes spin 
effects in the case of polarized tau decays. These simulations necessarily include form 
factors that model the resonant structure of intermediate hadronic states, such as p, u> 
or ai mesons. The combination of Monte Carlo simulation and experimental studies thus 
allows us to test the model input and leads to additional information about hadron physics 
at low energies. 

For a few channels only, and in a limited kinematical range, the form factors can be 
predicted from a firm theoretical basis. In many cases additional input, such as vector 
dominance and phenomenological parametrizations are required. At present only final 
states with up to four pions can be simulated on the basis of realistic form factors [3], 
which have their basis in chiral Lagrangians and vector-dominance models. In the present 
paper we describe an extension of TAUOLA to five-pion final states. The model for 
the amplitude is based on the observation that the 27r~7r + 27r° decay is dominated by 
the um~7T° channel, with only about 20% left for the remainder, which does not exhibit 
sharp resonance structures. The amplitude is therefore constructed to accommodate this 
feature. It includes the dominant W* — > u2n amplitude and an amplitude of the form 
W* — > ai(— > pvr)/ (^ 7T7r) with various charge assignments. The second amplitude is 
also used to describe the two remaining charge combinations, 7r _ 47r° and 3tt~2'k + . (The 
charge-conjugate combinations for the r + decay are considered in parallel.) In the narrow- 
width approximation, important technical tests of the generator can be performed. The 
details of the new amplitude (effectively the transition matrix element of the hadronic 
current, often just denoted "current") and of the phase-space generation are described 
in section 2. Tests of the program and results for some characteristic distributions are 
collected in section 3. Section 4 contains our summary and conclusions. 

2. Hadronic currents and phase-space generation 

For each decay mode the basic ingredients for TAUOLA can be grouped into three 
parts: (i) the phase-space generator for a fixed number of particles in the final state, (ii) 
the algorithm for the calculation of the spin-dependent matrix element from the hadronic 
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currents and from the properties of the electroweak interaction between the W boson, the 
r lepton and its neutrino, and (iii) the hadronic current per se. 

The structure of the program is decscribed in detail in refs. [1,3] (see also [4]) and 
there is no need to recall it here. The relation between the hadronic current and the 
matrix element remains the same as for the two-, three- and four-meson case presented 
before. 

Some comments are required concerning the five-pion phase-space generator, although 
its structure is very similar to the well-documented case of four pions. The number 
of generated angles and invariant masses is, of course, increased. Since no extensive 
studies of amplitudes with narrow resonances are foreseen for the moment, the phase- 
space presampler is only prepared to compensate for possible sharp peaks in the invariant 
mass of the five-pion system as a whole and in one subsystem of the three-pion (single 
channel) only. To improve the efficiency, this mode of generation is merged with the flat 
phase-space subgenerator, again according to the rules described in detail in [1,3]. 

The most interesting ingredient necessary for this generator is the matrix element of 
the hadronic current. For the decay into an n-pion state, it is given by 

■^(gi,92,-,5n) = (7r(?i),7r(?2),-,7r(?n)|^(0)|0), (1) 

where the same letter J has been used for the operator and its matrix element. The decay 
into an odd number of pions proceeds through the axial part of the current only. 

We now construct the amplitude as a sequence of (partly virtual) resonance decays and 
transitions, and concentrate, in a first step, on the most complicated channel 27r~7r + 27r°. 
From experiment we observe that the decay through u seems to constitute the dominant 
mode in this case, with a branching ratio Br(r — > h~Lun°) = (4.4 ± 0.5) x 1CT 3 , to be 
compared with Br(r -> h-h-h + 2n°(exl.K ,uj,t])) = (1.1 ± 0.4) x 10~ 3 [5]. The first 
part will be implemented in a first step (current A), and the remainder (current B) 
afterwards. Normalizations will be appropriately set, as can be seen in Section 3. From 
isospin conservation — and the fact that the total hadronic system has isospin one, and 
the u>, however, isospin zero - - it follows that the two pions are also in an isosp in- 
one configuration, corresponding to a state with p-meson-like quantum numbers. We 
thus have to find a convenient description for the amplitude that specifies the transition 
W*(Q) — > oop. Let us denote the polarization and momentum of the p meson by e p and 
p p , and similarly for the u meson. Three possible amplitudes describing the decay of a 
1 ++ to two 1 states can be constructed: 



Fp = (tp,£u),Pp PwjAO (cpjC^jPp Pu>,Q)~q^ 

Fp = (e p ,eu,Pp,Puj)((Pp-Puj)»- (Pp-P^Q^^), (2) 

each of which is, of course, multiplied by a function of Q 2 , where Q = £V q t . For simplicity 
we adopt amplitude F 1 , which depends on the lowest power of the relative momentum 
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Pp — Puj- Furthermore we multiply this amplitude by a Breit-Wigner factor 

simulating the a% enhancement. For the constant Co, which is determined by the product 
of W*a>i- and a i cup-coupling we adopt the value Co = 3, so as to reproduce the desired 
branching ratio of 0.4%. Furthermore we take m a = 1.26 GeV and T a = 0.4 GeV. With the 
symbol (e p , e w , Q, p) we denote the totally antisymmetric Levi-Civita symbol, contracted 
with three four- vectors e p , e w (polarization vectors for p and u), and Q; the last index p 
remains open. 

The amplitude for the "subsequent" decay of the virtual p is given by 

M^ p = g^(q 4 - q 5 Y, (4) 

which leads to the decay rate 



p ~ 48tt m 2 ' [) 



To reproduce the input parameters m p = 776 MeV and T p = 150 MeV the value g pwn = 6.0 
is adopted. For the u we adopt the corresponding decay chain u — > np, where all three 
Ti and p charges contribute with equal weight. We start from an u-p-n coupling of the 
form 

■M-UJpTT ^fuipwi^UJl PuJl €pl <?7r)? (6) 

where the coupling f upn has dimension (mass) -1 . Including the p decay according to 
eq. (Jl]), and taking the antisymmetric isopin wave function of the three points into account, 
which implies an antisymmetric momentum wave function, we finally arrive at 

Mua, = ^f^(e., gi, q 2 , q 3 )(BW p ( Sl ) + BW p {s 2 ) + BW p (s 3 )), (7) 
m p 

where qi, q 2 , q 3 denote the momenta of 7r + , 7r - , 7r°, and s\ = (q 2 + q 3 ) 2 etc. This leads to 
the decay rate 

x( Sl s 2 s 3 - ml(Q 2 - ml) 2 ). 

If the sum of the three p Breit-Wigners is dropped (i.e. replaced by a constant taken to 
be 1) and the pion masses are set to zero, one finds 

r - - - m 7 ^^P^Sptttt) 2 / q x 

u 3 128(2vr)3 m ^ 120m* ' W 
3 



This formula will be useful for cross checks of the program. Whenever numerical values 
are required, we take /^ p7r = 0.07 MeV -1 , that is only the first significant digit, g p7rn = 6 
was defined earlier, m u = 782 MeV and T w = 8.5 MeV. 

The full five-pion amplitude is obtained by including the p and to propagators 

</>i,feg 3 ,g4,g 5 ) = J-^f^BW a {Q 2 ) BW p {{q A + q 5 ) 2 )BW u (( qi + q 2 + qs?) (10) 

0, 04 - Ob, a, Q)(a, qi, q 2 , q 3 )(BW p (s 1 ) + BW p (s 2 ) + BW p (s 3 )). 
The pictorial illustration of this decay amplitude is shown in Fig. la. 




Figure 1: Dominant decay amplitude for the decay of r into five pions through an uj plus 
a p resonance (a) and through an f plus a\{— > pn) (b). 

This formula is directly applicable to the case of a narrow u, if we identify, as stated 
above, qi, q 2 , q 3 with the momenta of 7r + , 7t~ and 7r° from the u decay, q^ and q$ with the 
momenta of the remaining ir~ and 7r°. Without this approximation the symmetrization 
between q 2 and q± on the one hand, and q 3 and q§ on the other hand has to be performed 

J p ' t0ta \qu 02, 03, 04, Ob) = J p (0l, 02, 03, 04, 05) + j£(0l, 04, 03, 02, 05) + 

J p (0l, 02, 05, 04, 0s) + tfiOl, 04, 05, 02, 0s) (11) 

and the phase space has to be divided by 4 to take the identical particle statistical factor 
into account. 

For a test of the proper implementation of the matrix element (eq. ()10J0. the narrow- 
width approximation for oj, p and uj is employed. (For this test (m p + m w ) < m a must be 
assumed.) In this approximation the rate derived from the current J A can be integrated 
analytically. (Again the sum of three p Breit-Wigner amplitudes in the uj decay is replaced 
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by 1 and m n is set to zero for this test.) 

?A ?e = 7—^(1 a 1 ) ( 1 + 2 ^ L ) c o-p — f(™ ai ,m p ,mJR 

4 m~ V / V m; / 1 ai ^ 

2p / p 2 p 2 



f( 2 2 2^ _ n Z P ( P , P "o\ 



p 

A 1 / 2 (m^ ,m 2 m 2 ,) 

p = V ai ' p; ^ , (12) 

2m ai 

The factor R = tt-^t- consists of the product of the partial widths for p and uj as given 
by eqs. (0) and 0, divided by total widths used as numerical inputs in the Breit-Wigner 
amplitudes, and the decay rate is normalized to T e = T(r — ► eu e i/ T ). 

For a test of the generator, in particular of the phase-space integration, the following, 
totally unphysical form of the current 

J% — c i Qn 
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Cl = -^-4!(4tt)V20 (13) 



was used, with the pion mass again set to zero 1 . The analytical result, 

T c /T e = cos 2 9 Cabbibo = 0.950625, (14) 

is well reproduced by the generator. The numerical results of the second test, based on 
the narrow- width approximation, of eq. (J12)) . will be discussed below. 

The decay mode into ujtctt is obviously only possible for the 27r°27r~7r + final state. In 
contrast the amplitude introduced in the following will contribute to all three charge com- 
binations 7T~47r°, 27T~27r°7r + and 37r~27r + . For this second amplitude we use the transition 
of the virtual W into an a\, with subsequent transition into a\ plus two pions in an isospin 
zero, angular-momentum zero configuration, parametrized by the broad fo resonance. For 
the a\a\fo coupling the simplest Lorentz structure with minimal momentum dependence 
has been adopted. The ai subsequently decays into a p meson plus a pion, with equal 
amplitude for the two charge modes n~2n and tc + 2tc~ , similar to the parametrization of 
the tau decay into three pions. 

This leads to the following current: 

Ju(Pl,P2,P3,P4,Ps) = , 2 o faafff™9apir9pKK X 

BW a (( Pl +P2+P3+P4 + P5) 2 )BW a (( Pl +p 2 + p 3 ) 2 )BW h ((p 4 + p b ) 2 ) x 

" 9,u) ( / 2 | Pl 7 3 \ 2 fa+^+P3r " (Pi -Psr)BW p (( Pl + P3 ) 2 ) 

V Q 2 J \ {Pl+P2+P3> J 

+ (1 



(15) 



1 As an alternative, in principle we could also use C2 = 4!(47r) 3 /Q 3 and obtain the same result. In 
such a case, numerical-stability problems could appear as, for masless pions, may approach integrable 
infinity within the allowed phase space. 
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with the momentum assignments d\(— > 7r-(pi)vr_(p 2 )7r+(P3)) + f(—> 7r o(P4) 7r o(P5))- The 
pictorial illustration is shown in Fig. lb. (In the program, we have used: rrif = 0.8 
GeV, = 0.6 GeV, G apw = 6, f aa f = 4, ff nn = 5 and c = 4.) The last constant was 
introduced to normalize the branching ratio for this channel to 0.11%. The amplitude is, 
by construction, symmetric under the exchange p\ vs. p 2 and p^ vs. ps, as requested from 
Bose symmetry for the fo and a\ decays. Alternatively, the two 7r° may originate from 
the a\ with the fo then decaying into 7r - 7r + . In this case the symmetrization with respect 
to the momenta of the ir~ has to be performed explicitly. For consistency we have to 
adopt the same momentum assignment as before: 7r + (q , i)7r~(g2)7r (5 , 3)7r~(g4)7r (g 5 ). The 
properly symmetrized amplitude thus reads: 

jZ> tatal (q 1: q 2) q 3: q h q 5 ) = J B (q 2} g 4 , q 1: q 3} q 5 ) + (? 3 , 45 , 42, 4i, 44) + ^(ffa, ?s, ?4, 9i, 

(16) 

We include the same statistical factor of \ as before into the normalization of the phase 
space. The relative rate [f(-\ — ) + a(00)] : [/(00) + a( — h)] = 2 : 1 is recovered in the 
narrow width approximation for a\ and fo- At present there is no analytical benchmark 
for the overall normalization available for this channel. 

The extension of this model to the description of the remaining charge configurations 
is straightforward. Let us start with tt~4tt° and adopt the following momentum assign- 
ment: vr~(g 3 )7r (gi)7r (g2)7r (94)7r (Q , 5)- Using the definition (fTB} for J B as in eq. (JIEJ), the 
amplitude has to be symmetrized with respect to the momenta of 7r°. Since J B is already 
symmetric with respect to the first two and the last two momenta, only 6 out of the 4! 
permutations have to be considered. This leads to the current 

Jj°~ 00 (gi, 42,43,44,45) = ^(41,^2,43,44,45) + ^(45,42,43,44,40 + ^(42,44,43,41,45) 

+ Jj? (4i, 44, 43, 42, 45) + ^(4i, 45, 43, 44, 42) + ^(44, 45, 43, 4i, 42)- 

(17) 

For the evaluation of the rate the statistical factor of 1/4! must be included. 
In complete analogy we obtain for the 37r - 27r° mode 

^ _++_ (4i, 42, 43, 44, 4s) = ^(4i, 42, 43, 44, 4s) + (45, 42, 43, 44, 4i) + J*(qi, 4s, 4s, 44, 42) 

+ J* (4i, 42, 44, 4s, 4s) + (45, 42, 44, 4s, 4i) + (4i, 45, 44, 4s, 42), 

(18) 

where the momentum assignment vr _ (g 1 )7r _ (g 2 )7r + (g3)7r + (g4)7r _ (g 5 ) has been adopted and 
a statistical factor l/2!3! = 1/12 is used in the evaluation of the rate. 

For the "non-a;" decays and in the narrow-width approximation, this ansatz predicts 
the following abundances of the subchannels: 

/(00)a(- - +) : /(-+)a(00-) : /(OO)a(OO-) : /(-+)a(- - +) = 1 : 2 : 1 : 2 (19) 

and the following relative rates 

27T°27r"7r + : tt^tt : 3^2^+ = 3:1:2. (20) 
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Note that this prediction is specific to the resonance and isospin structure of the model. 
Once more experimental information on the five-pion channel will be available, more 
elaborate possibilities can be considered. 

3. Results from the Monte Carlo Program 

The new version of TAUOLA includes new channels, numbered 31 to 35, which will 
be discussed in turn. Channels 31, 32 and 33 refer to 27r°27r - 7r + , channel 34 to tt^Att 
and channel 35 to 37r~27r + . Channels 31 and 32 are implemented for tests, channels 33, 
34 and 35 for physics simulations. 

Let us start with channel 31, which is based on current J A , without symmetrization, 

as defined in eq. (JTUJ). (The order of the pions generated by TAUOLA is ( h 0). This 

ordering is necessary for the test and the requirements of the presampler running in the 
narrow- resonance mode, which is only partly optimized. The result, based on the default 
parameter values listed after eq. Q, is given in the first line of Table 1. This channel 
is then used to test the program against eq. (fT2"|) . which is valid in the narrow-width 
approximation for p, u, and a, using massless pions. This simulation was quite demanding 
on the phase-space generator. A precision of about 1% only was reached. Because of the 
complex structure of the resonances a%, p and u it was impossible to generate events 
efficiently. The p Breit-Wigner was only generated from a flat distribution. The variance 
of the generated raw sample was deteriorating quite fast with decreasing width, not only 
rendering the generator slow, but also risking to arrive at numerically unstable results. 
We nonetheless completed this test for several choices of masses and widths of a±, p and 
to. A typical result is listed in Table 1, line 4, where = T a = 1 MeV, T p = 5 MeV, 
nip = 373 MeV, was adopted for the parameters, and other parameters were left at their 
default values, in particular m w = 782 MeV and m a = 1260 MeV was kept. The remaining 
2.5% difference between Monte Carlo and analytical calulation is due to contributions 
from the tails, which remain large, even for the extremely narrow widths adopted in this 
example. This illustrates the importance of non-resonant contributions and interferences, 
in particular for realistic values of the parameters. 

Current J c , which also serves to test the generator, is implemented in channel 32. In 
the case of massless pions, eq. (|T4]) is reproduced with a precision better than 0.1% (see 
line 5 of Table 1). This result is stable and independent of the choice of options for our 
phase-space presampler, providing an important test of its correct operation. The result 
for J c , taking the pion mass into account, is listed in line 3. 

In the second line of Table 1 the result is listed for current J B , as defined in eq. (fT3j) , 
again without symmetrization. 

The results for physical input parameters, and with properly symmetrized amplitudes, 
are listed in Table 2, column 4. The first three lines refer again to 27r°27r~7r + , as imple- 
mented in channel 33, and are based on eqs. (fTTj) and ()16|) . For the order of the pion 
momenta in the program the convention 2 (H 00) has been used. The first line is 

2 Note the difference in order with respect to Section 2 and channel 31. 
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Current 


T x /T e (TAUOLA) 


Yx/Y e (analytical) 


Comment 


A 
B 

C 


r\ r\r\ A r\ A i o — 

0.02404 ±3 xlO 
0.00928 ±9 xl0~ 6 
0.10690 ±2 xlO" 4 




without symm. 
without symm. 
massive pions 


A 
C 


1.38xl0 7 ± 1 xlO 5 
0.95030 ± 3 xlO" 5 


1.41 xlO 7 
0.950625 


massless pions 
massless pions 



Table 1: Test results of the generator for unphysical choices of parameters; see text. 

based on the ujtitt channel ( j A ; total ) alone, and the symmetrization evidently increases the 
result by about +5%. The second line is based on j B ' total alone, and the symmetrization 
decreases the result by 11%. Adding coherently j A ' total and j B < total one arrives at the 
result given in line 3. Keeping in mind our previous experiences with interferences, it is 
amazing that the result is exactly the sum of the two previous entries. This implies that 
the two amplitudes do not lead to interferences in the rate in any significant manner! But 
of course, perfect agreement for all significant digits is due to statistical fluctuations. 

The charge configurations tt~Att° and 37r~27r + (with the momentum ordered (—0 0) 

and ( — h H ) in the TAUOLA output, again differently as in Section 2) correspond 

to channels 34 and 35 respectively. The results are based on the amplitudes given in eqs. 
(117)) and (J18j) and the numerical predictions are listed in lines 4 and 5. These numbers 
are contrasted with the experimental results listed in Table 2, column 5. The result for 
27T~7r + 27r° as obtained with J A ' total alone is compared with the mode h~um°, the result for 
j_B ; totai alone with the mode h~h^h + 7i 0, K (exl. K°,cu,r]). The 7r _ 47r° mode is compared 
with h~4n (ex K°,r]) and 3n~2iT + with h~h~h~h + h + . Let us emphasize again that the 
rates for the two submodes of channel 33 are fitted to the data, the rates for channels 34 
and 35 are then the result of our ansatz. 

The relative rates as displayed in lines 2, 5 and 6, column 4, are seemingly in contradic- 
tion with the expectations based on the narrow- width approximation for the intermediate 
states, eq. (|20|) . We therefore perform the simulation also in the narrow- width approxima- 
tion, with the following parameters: m/ = 0.45 GeV, Yf = 10 MeV. For the a saturating 
intermediate 3-pion state, we took m a = 1.16 GeV and Y a = 10 MeV, while width and 
mass of the a± coupled directly to the r lepton were kept at their standard value; this 
was also true for the mass and width of the p meson. Indeed, in this completely fictitious 
case, eq. (120)1 is recovered within the statistical error (column 6). 

Eq. (120)) is only reproduced for very small widths of the resonances. For realistic 
parameters, the predictions for ratios of branching fractions are significantly affected by 
interferences, and this applies even more so for differential distributions. This observation 
has to be taken into account from the very beginning in the formulation of a realistic 
current and the construction of Monte Carlo programs where cascade decay chains must 
be implemented. 
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Channel 


Final state 


Current 


T x /T e x 10 3 


T x /T e x 10 3 


Tx/r e 










(TAUOLA) 


(experiment [5] ) 


(narrow width) 


33 


27r-7T + 27r° 


A 


total 


25.30 ±0.1% 


25 ± 3 




33 


27T-7r + 27r° 


B 


total 


6.05 ±0.2% 


6.2 ±2 


2192 ±22 


33 


27r-7r + 27T° 


A+B 


total 


31.35 ±0.1% 


31 ±2 




34 


Ti 'Air 


B 


total 


9.37 ±0.1% 




788 ±4 


35 


37T-27T+ 


B 


total 


11.03 ±0.1% 


4.6 ± 0.3 


1469 ± 7 



Table 2: Test results of the generator for realistic choices of parameters; see text. 

4. Summary 

In this note, recent developments in TAUOLA for r decays into five mesons are described. 
For the 2n^n + 2n° mode, two amplitudes were introduced. The first is motivated by the 
experimental observation of a dominant ujtttt channel, the second corresponds to a non-u; 
contribution and is parametrized by the transition through an intermediate a\ + fo com- 
bination. This second amplitude was also used to predict the two remaining charge com- 
binations for r decays into five pions. Both amplitudes are implemented into TAUOLA. 
A third current, J c , was introduced for tests of the phase-space generator. With the help 
of these currents basic technical tests of the algorithm were successfully completed. No 
sophisticated fits to the data were attempted. The model amplitude can, however, be 
used as a starting point for more detailed theoretical and experimental investigations. 

We observed that effects due to interferences between different resonant amplitudes 
are often large. This is due to various threshold effects and may be specific for the choice 
of our amplitude. It may result, for example, from those various threshold effects or from 
simplifying assumptions adopted for our currents, such as the choice of constant instead of 
running, final-state mass-dependent widths. For different parametrizations, the influence 
of interferences on the total rates could be significantly smaller. 

This isssue is of great practical importance for the design of a generator used in fits 
to data during their analysis, when the general form for the currents often has to be 
changed several times. If the constraint of zero interference between various currents, 
standing for different subresonances, was imposed on the algorithm, as was the case e.g. 
with EURODEC [8,9], the fitting procedures would become unnecessarily difficult (even 
impossible), or might even not lead to a correct description of the amplitude. Therefore, 
distinct currents for each final state, which are labelled by stable decay products, will be 
also used in the future. Resonances such as p, a±, u are can only be used as intermediate 
states in the amplitudes. 

Finally, let us point out that the extensions for TAUOLA as presented in this pa- 
per are available from the directory tauola-BBB of the standard distribution system for 
TAUOLA versions. The usage of the directory tauola-BBB is not documented yet; it is, 
however, identical to the one of the tauola-F subdirectory as explained in [4]. The most 
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recent version of the code is available from the web page [7]. Test results for differential 
distribution, obtained with the help of MC-TESTER [10], are available from the web 
page [11]. 
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